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Abstract 

The electronic structure of lanthanide and actinide compounds is often characterized by orbital 
ordering of localized /-electrons. Density-functional theory (DFT) studies of such systems using the 
currently available LDA+?7 method are plagued by significant orbital-dependent self-interaction, 
leading to erroneous orbital ground states. An alternative scheme that modifies the exchange, not 
Hartree, energy is proposed as a remedy. We show that our LDA+^J approach reproduces the 
expected degeneracy of f 1 and certain f 2 states in free ions and the correct ground states in solid 
PrC>2. We expect our method to be useful in studying electronic excitations and entropies in /- 
and he&vy-d elements. 

PACS numbers: 71.15.Mb, 71.27.+a, 71.20.Eh 



I. INTRODUCTION 



Interesting physical phenomena associated with the strongly correlated /-electrons in 
lanthanide and actinide compounds continue to attract lively interest!-^. Strong on-site in- 
teractions between the /-electrons in these materials present serious challenges to modern 
density-functional theory (DFT) based electronic-structure techniques, causing most approx- 
imate functionals, such as the local density (LDA) or generalized gradient approximation 
(GGA), to fail qualitatively. To overcome the deficiencies of the LDA/GGA in studying 
/-element compounds, several recent studies have employed the self-interaction-corrected 
LDA 3 e.g. in Refs. 11151151 the hybrid functional method 33 in Refs. IMTUlfTTllT2lfT31 or the dy- 
namical mean- field theory (DMFTjPin Refs.QSCSI The LDA+U method 13 has emerged as 
a well-established model to deal with strong electron correlations in d- and /-systems, com- 
bining high efficiency with an explicit treatment of correlation within a Hubbard-like model 
for the localized electrons. This method has been very successful in transition metal oxides 
(for a review see Ref. ITSl) and has yielded promising results for band gaps in / systems ^ 3 * 19 * 20 !. 
However, systematic studies of its effectiveness remain inconclusive, with issues of orbital 
ordering^ and multiple self-consistent solutions attracting heightened attentio n 119 1 22 1 23 1 24 12 ^ 1 . 

Here, we show that the currently popular versions of LDA+f/, by Liechtenstein and co- 
workers^ 3 and by Dudarev and co-workers^, respectively, encounter serious difficulties in 
/ systems due to large orbital-dependent self-interaction (SI) effects, which result in an 
unphysical splitting of up to 0.4 eV between degenerate f 1 multiplets. Since the SI errors 
(SIE) are typically larger than the crystal field (CF) splitting energies, and comparable to 
the strength of the spin-orbit coupling (SOC), they lead to qualitatively incorrect electronic 
ground states in solids. We propose a new, orbital SI free form of the LDA+[7 method 
that leaves the LDA Hartree term intact and only replaces the LDA exchange with the 
Hartree-Fock exchange. In our method, the Hartree-Fock exchange term cancels the LDA 
self-interaction energy to a high degree of accuracy, ensuring near-degeneracy of real- and 
complex-valued orbitals in free ions and correctly reproducing the Tg ground state and 
Tg — > Tj excitation energies in the Pr02 solid. The accuracy of this functional is sufficient 
for evaluating high-temperature electronic entropies of / electron systems. 
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II. METHOD AND COMPUTATIONAL DETAILS 



All DFT calculations were carried out using the VASP package^^ with projected aug- 
mented wave (PAW) potentials 33 , energy cutoff of 450 eV, and without any constraint sym- 
metry or ionic relaxation. For free ions, a 12 Acubic cell containing one ion and uniform 
compensating background charge were used. For the PrC>2 solid, we consider a primitive cell 
of the fee supercell (lattice constant of 5.386 A^B). The term "LDA+t/" is used irrespective 
of the xc functional since the LDA and GGA results are found similar. Each calculation was 
initialized in a specific atomic orbital and self-consistently converged to either states very 
close to the initial orbital with the results reported, or distinctly different states with lower 
energy. SOC was excluded from the calculations unless its inclusion is stated explicitly to 
make realistic comparison with experiment. Finally, we fix the U parameter in the LDA+[7 
method to 6 eV and leave the discussions of this choice to the end. 

A. Aspherical self-interaction error of LDA+[7 for /-electrons 

We begin by showing that the conventional LDA+?7 approach fails to reproduce the 
degeneracy of different \m) orbitals of f 1 ions. First consider real orbitals with angular 
dependence of real yf m = ^/2?RY 3m without spin-orbit effects to simplify the presentation of 
our method. Complex orbitals and SOC are discussed later. Fig. [T] shows the energies of 
different y^ m orbitals (with the exception of y^, which converges to y^ 2 ) in several lanthanide 
and actinide ions calculated using the LDA+£7 scheme of Liechtenstein et a/P^ with J = 0.5 
and U = 6 eV. Contrary to the expected degeneracy, the energies of the different y§ m orbitals 
differ substantially, up to 0.4 eV, and the y§± orbital was found unstable and converged to 
7/32- Varying J between (i.e. the Dudarev scheme^) ~ 1 eV changes the results by only a 
few meV. 

The above results demonstrate that the conventional LDA+£7 approach commits large 
errors of up to 0.4 eV/electron in the predicted relative orbital energies of / electrons. 
To understand the reasons for the unphysical splitting of the f l states, we examine the 
conventional LDA+£7 total energy functional 

E LDA+U = E LDA + Eu _ ^ (1) 

where the LDA description of the on-site interaction, approximately represented by the so- 
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called double-counting term E^c, is replaced with a Hubbard-like E\j. The latter is essentially 
the Hartree-Fock energy, expressed in a rotationally invariant form by Liechtenstein et a/P^ 
as a sum of the Hartree (H) and exchange (X) terms, Ey = En + E x , where 

Er = ^^(m,m"\V ee \m\m'")n mm m m >> m m, (2) 

{m} 

E x = ~ ^ (m,m''\V cc \m''',m')n^ nmf n^ n ,, ml , l . (3) 

{m},cr 

The on-site density matrix n a mml is obtained by projecting the Kohn-Sham orbitals ip^ of 
occupancy onto atomic states \nlm(m')) 

<m> = Y,f«(V>lrn'){nlm\r a ), (4) 

a 

while the Slater integrals {mm > \V ee \m"m" > ) are evaluated in terms of the Gaunt coefficients 
and the screened Coulomb U and exchange J parameters (the diagonal m — m' — m" = wl" 
terms are given in Table |l]). A simplified version by Dudarev et alW^ adopts the J = limit. 
States with real n°, are referred to as "real". 




FIG. 1: LDA+C7 total energies for different orbital filling of f l ions with the Liechtenstein scheme^ 
relative to y^ 2 . 

For a free f 1 ion, the Hartree-Fock energy Ejj in Eq. [T] naturally vanishes, while E^c in 
the Liechtenstein and Dudarev schemes depends only on the number of electrons, N u = 
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TABLE I: Hartree energy En, Eq. [2] and LSD exchange energy -E^ SD , Eq. ([7]), for one /-electron 
in orbitals with real (y/^,^o) an d complex (YJ m for m > 0) angular wavefunctions. 

X] m ^mm' anc l n °t 011 t ne type of the occupied orbital. Therefore, Eq. [I] becomes 

E LDA+U = E LDA + congt _ ^ + congt _ 

In the above approximation we assumed 1) the LDA exchange is not sensitive to orbital 
filling (more on this later) and 2) the Hartree energy difference comes mainly from the 
on-site Hartree term i?H of eq. [2] The resulting error in the relative orbital energies is 
then entirely due to the orbital-dependence of the SIE of the LDA, which is reflected in 
£"h. To see the validity of our argument, we list in Table |T] the on-site En calculated from 
eq. [2] for atomic orbitals; these expressions are expected to closely approximate the SI for 
localized orbitals in the LDA+U. Even though is identical for all real p or d orbitals, 
it is orbital-dependent for / multiplets, and in all cases splits the SI energies of real vs. 
complex orbitals. The predicted ordering of En is 2/32 < 2/30 < Z/31 < 2/33, in agreement with 
the LDA+U results shown in Fig. [TJ demonstrating that the unphysical splitting of f 1 states 
in conventional LDA+f/ is due to orbital-dependent SIE. Note that with real orbitals the 
problem of orbital-dependent SIE does not affect p or d electrons. We will show later that 
complex p and d orbitals are affected. According to Table [TJ the SIE is proportional to J; for 
typical values of J in the range of 0.1 to 1 eV, it is comparable to or even larger than other 
important on-site effects, such as CF and SOC, which can lead to qualitatively incorrect 
predictions of electronic ground states in solids by the current LDA+f/ methods. These 
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deficiencies of the conventional LDA+U approach can be traced back to its treatment of the 
Hartree and exchange energies. The LDA+U approach replaces the LDA Hartree energy 
with an on-site model expression given by Eq. [2j Even though the E^ term is capable 
of reproducing the correct orbital energetics, the LDA+U double-counting energy E dc is 
orbital-indepedent and fails to properly account for the orbital-dependence of the LDA SIE 
in open-shell systems. Similar considerations hold for the orbital-dependence of the LDA 
exchange energy, which is mainly sensitive to the choice of real vs. complex orbitals (see 
Table [I]); this factor acquires importance in systems with strong SOC, when the orbitals 
with a definite value of the total angular momentum J are necessarily complex. 



B. Reformulated LDA+U 

To correct the orbital-dependent SIE in the Hartree and exchange terms, we propose a 
new formulation of the LDA+U method by modifying only the exchange term of the LDA: 

E LDA + U = E LBA + Ex _ ^ (5) 

where the orbital-dependent Hartree-Fock exchange Ex of Eq. ^ contains a term that 
approximately cancels the SIE in the LDA Hartree energy; the remainder of the LDA Hartree 
energy is exact by definition and therefore left unmodified in our approach. The exchange 
double-counting term E^cx accounts for the LDA exchange energy and is given by a linear 
combination of the exchange double-counting in the Liechtenstein scheme and the on-site 
local-spin-density (LSD) exchange: 



E dcX = $^iV CT + JN°(N° - 1)] + cE™ D , (6) 

cr 
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where c is the interpolation coefficient, p a is the charge density of spin component a, which 
can be obtained from the on-site occupation matrix n a mml as well as radial function Ri(r) and 
spherical Y/ m (f2), K is the LSD exchange strength parameter, and p represents the angular 
part of p. Only the Elf D term in Eq. (JsJ) is orbital-dependent. The linear interpolation is 
conceptually similar to hybrid functional approaches and serves the purpose of subtracting 
the orbital-dependence of the LDA exchange energy. The potential corresponding to the 
correction energy Ex — E dcX , obtained by differentiating with respect to the on-site density 



matrix n mm >, is then 



2c / An ^ 1/3 



3 \2l 



-j k dnip°(n)} 1/3 Y lm Y lm , 



+ (1 - c)(^y^ + rfj)8 mm , - (rn, m"\V ee \m"'m')n a m „ m „, (8) 

It is possible to reduce the number of independent parameters by requiring that E x — E dcX 
vanishes for full /-shells (n^ m , = n l mm , = 5 mm >), 

Ex - £dcx = -c(2l + 1)(U + 2U) + c(2l + l)K = 0, 

which gives 

K = U + 2U. (9) 

The main advantage of Eqs. ([5])- ^ is that the LDA self-interaction energy is canceled by 
the corresponding exchange term in Ex- As a result, the proposed method is self-interaction 
free to high accuracy. 

III. RESULTS AND DISCUSSIONS 

In this section, we analyze the parameter dependence of the proposed method and then 
presents results for the example of PrC>2 solid. 

A. Determination of parameters to remove aspherical SIE 

We demonstrate orbital degeneracy for free Pr ions with one and two /-electrons. Fig- 
ure |2^i displays the energy of Pr 4+ in real atomic orbitals calculated with our method (as- 
suming c = 0) as a function of the exchange parameter J. At J = 0, a splitting of more 




J/eV 

FIG. 2: LDA+C/ energy of the Pr 3+ and Pr 4+ ions as a function of J for c = calculated with 
our method, a) f 1 in orbitals y| m , with the optimal J region magnified in the inset; b) f 2 in three 
degenerate (in terms of Ejj) two-electron states. 

than 0.3 eV is found, similar to the behavior of the original LDA+£/ in Fig. [TJ The splitting 
is reduced by increasing J and at the optimal value of J° = 0.783 eV, it is less than 40 meV, 
i.e., the four real orbitals y§ m are almost degenerate. The y 3l orbital can only be stabilized 
in the vicinity of J°, relaxing otherwise to the more stable y 32 or y^ 3 . Hence, just one point 
for y§y is shown in the inset of Fig. |2^i. 

The energy of the Pr 3+ ion (f 2 ) is shown in Fig. [2]b (also at c = 0). Consider three distinct 
f 2 states with S = 1 and degenerate Hartree-Fock energy E v . Using the basis denned by 
real-valued spherical harmonics, {y{\ m \ = V2$sYi\ m \(— I < m < 0), Yi , y^(0 < m < /)} 
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0.0 0.2 0.4 0.6 0.8 1.0 

c 

FIG. 3: LDA+J7 energy of the Pr 4+ ion as a function of c at fixed J° = 0.783 eV calculated with 
our scheme, including both real- and complex-valued (no SOC) orbitals. 

(shown for I = 3 in Fig. the first of these states, designated by 0i 3 , has electrons in 
orbitals y 31 and y| 3 , or n^ m , = except n| x = nl 3 = 1, while the other two f 2 states, 
designated by 0i 4 and 0i 5 , correspond to n u = n\ A = 1 and n| x = n' 55 = 1, respectively. 
Their angular wavefunctions are shown in Fig. ^p. Similar to the f 1 case, the energy splitting 
is large at J = and gets reduced to less than 30 meV at the optimal value J°. Note that 
015 can be stabilized only for J > J°. 

So far, we have used c = 0, assuming that the LSD exchange functional is insensitive 
to the orbital and can be ignored. The lower part of Table [I] proves this assumption for 
the real orbitals: E^f D varies by less than 0.02K. However, Table [i] also shows that 
of complex orbitals is substantially lower (by ~ 0.37^), indicating a large lowering of the 
exchange energy in states with nonzero orbital current. Since ~ — p 4 / 3 is concave, it 

favors inhomogeneous charge distributions (such as real orbitals compared to complex ones) 
and therefore the LDA exchange energies in Table [I] of real orbitals are lower than those for 
complex orbitals. The difference may play an important role in systems with strong SOC, 
when the resulting electronic states are complex combinations of real yimS with the orbital 
angular momentum unsuppressed. In Fig. |3j we show the dependence of the energies of 
real- and complex- valued orbitals for Pr 4+ on the mixing coefficient c in Eq. (TH), using the 
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optimal value of the exchange parameter, J°. It is seen that at c = 0, the energies of real 
and complex orbitals differ by more than 0.2 eV due to their different LSD exchange, and 
the spurious splitting is minimized to approximately 70 meV at the optimal c ~ 0.6. 

In our approach, the J and c parameters are a priori determined by the physical require- 
ment of degeneracy once the U parameter is given (6 eV in this work). They hardly change 
when U = 4 eV is used, suggesting that our method is relatively insensitive to the choice of 
U. 

B. Eigenstates of Pr02 without SOC 
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TABLE II: Comparison of the LDA+t7 energy of PrC-2 and crystal field eigenvalues. Ground state 
energy is given in bold. 

Finally, we demonstrate the advantages of our method for extended solids by considering 
Pr02 in the cubic fluorite structure. The Pr 4+ ion is coordinated by eight oxygen atoms 
in a cube. Figure [4] shows the f 1 energy level splitting scheme in the presence of cubic 
CF and SOC. Without SOC, the cubic CF splits the f 1 states into the t\ u ground state 
and t 2u , d2u excited states (see Fig. |4^l,c) . Table [LT| lists the CF eigenvalues of these states 
(small 6th-order CF ignored), and the calculated LDA+U energies using the conventional 
approaches and our new scheme at the optimal values of J = 0.783 eV and c = 0.6. The 
conventional schemes predict orbital enegies that deviate dramatically from the expected 
CF order: the Liechtenstein approach predicts almost degenerate t\ u and t 2u , while t 2u is 
the ground state in the Dudarev method. In contrast, our new method successfully finds 
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the correct t\ u ground state. 



C. Eigenstates of Pr0 2 with SOC 
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FIG. 4: a) Schematics of f 1 energy levels (with multiplicity) split by SOC and cubic CF. The 
angular wavefunctions are shown for b) real- valued atomic orbitals ^3, CF eigenstates c) without 
and d) with SOC. 

The physics of orbital ordering in / systems is affected by strong relativistic effects^, 
necessitating the inclusion of SOC to make direct comparisons with experiment. Including 
SOC, our method predicts that the energies of the CF-degenerate Tg and r|, and the excited 
T 7 states in Pr02 (Fig. [4jd) are (reference), 69 and 142 meV, respectively. The spurious 
69 meV splitting between the two degenerate T 8 states is consistent with the accuracy 
shown in Fig. |3j Neglecting Jan- Teller lattice distortions and magnetic ordering effects, we 
estimate that the r 7 /T 8 CF splitting is between 73 and 142 meV, in good agreement with 
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the measured value of 131 meV from neutron diffraction^. 



D. Aspherical SIE in other methods 
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FIG. 5: Total energy of the Pr + ion for different filling of real orbitals as a function of cjexX) 
the fraction of exact exchange, with the hybrid functional (HSE06) method. Most calculations for 
complex orbitals Y^ m converged to very different states and are not shown. 

Our method bears some likeness to the hybrid functional approach. The difference in 
the latter is that the exchange interactions are calculated directly from the wavefunctions, 
with the amount of exact or Fock exchange (U/2 + a J for one localized electron in terms of 
LDA+C/) as well the replaced LDA/GGA exchange controlled by a fixed parameter czexx- 
However, oexx in the hybrid functional method is often system-dependent and fitted to 
experimental data, just like U in LDA+C/. For instance, Ref. [33] found that in /-ystems 
good results were obtained using 40 — 70% Fock exchange, while <i-systems typically require 
20 — 50%F^. However, such an ciexx m &y not necessarily lead to accurate removal of the 
aspherical SIE. Fig. [H] shows the energy of Pr 4+ ion as a function of ciexx calculated with 
the hybrid functional (HSE06pH Nearest degeneracy is obtained at ciexx ~ 85%. Given 
the sensitive orbital dependence of SI demonstrated in this work, in general the accuracy 
of hybrid functional calculations for /-electron systems may still suffer from incomplete 
removal of aspherical SIE. After the first submission of this manuscript, we became aware 
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that the idea of removing on-site E-r from LDA+?7 was previously proposed from a different 
perspective in Ref. EH1 in which the correction energy is independent of the orbital filling, 
an important different from our approach. Therefore, the method of Ref. ESis not expected 
to give accurate removal of the orbital-dependent SIE. 

IV. SUMMARY 

In summary, we have identified a serious problem in applying the LDA+[7 method to 
/-electron systems: the degeneracy of atomic orbitals is lifted, resulting in qualitatively in- 
correct electronic ground states and orbital excitation spectra. Aspherical orbital-dependent 
self-interaction is identified as the main source of error. To correct it, a new LDA+[7 scheme 
is proposed, which leaves the Hartree intact and only replaces the LDA exchange with the 
Hartree-Fock exchange. Our method has one adjustable parameter U, with the other two 
(J and c) being determined from the condition of orbital degeneracy in free ions. The 
computational expense is approximately the same as in the conventional LDA+C/, and very 
competitive compared to hybrid functional approaches^. We expect that our method will 
scale to large systems and will significantly improve the accuracy of first-principles studies 
of /- as well as heavy <i-systems with significant relativistic effects. Additionally, more ad- 
vanced methods such as GW and DMFT could benefit from the correct input ground state 
orbitals generated by our method. 
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